Mixed-Effects Models

Repeated Measures Model
Dataset: Mohlman, Jessica L. et al. (2019), Data from: Nonconsumptive effects of hunting on a nontarget game bird, Dryad Dataset
Northern Bobwhite quail

Disturbance through human hunting activity can significantly impact prey species through both consumptive and nonconsumptive effects. The nonconsumptive effects of rabbit hunting on Northern Bobwhite (Colinus virginianus; hereafter, bobwhite) could cause an increased perceived risk of predation by bobwhite during rabbit hunting events may elicit anti-predator responses, such as reduced movement away from the safety of cover.

Northern Bobwhite quail
bobwhite <- read.csv('bobwhite3.csv')
bobwhite$ID <- as.factor(bobwhite$ID) #make ID a factor
p1<- ggplot(bobwhite, aes(x=HuntDay, y=HW_Dist, group=ID, color=ID, shape=ID)) + 
  geom_point(size=4, alpha=0.6, position = position_dodge2(width=.33, preserve = "total")) +
  scale_y_continuous() +
  #geom_line() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title="Risk Behavior in Bobwhite During Hunting Season", x= "Hunting Season Species", y = "Distance from Hardwood Forest Cover (meters)")+ 
  theme_bw()+
  scale_color_brewer(palette = "BrBG")
p1
## `geom_smooth()` using formula 'y ~ x'

# ID: Unique ID given to each bobwhite covey tracked in chronological order.
# HuntDay: Denotes if it was a “Rabbit” or “Quail” (bobwhite) scheduled hunt day. 
# HW_Dist: Distance in meters a bobwhite covey was from hardwood habitat.
mixed_bobwhite <- lmer(HW_Dist~(HuntDay*ID)+(1|ID), data = bobwhite)
## Warning in as_lmerModLT(model, devfun): Model may not have converged with 1
## eigenvalue close to zero: 5.5e-09
anova(mixed_bobwhite) 
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## HuntDay    19204.9 19204.9     1   188 11.0107 0.001088 **
## ID          3168.5   633.7     5   188  0.3633 0.873151   
## HuntDay:ID 23341.2  4668.2     5   188  2.6764 0.023103 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Looking at all of the possible interactive effects, we can see that the only statistically significant effect is the interaction of hunting season day and how close the covey stayed in the shelter of hardwood forest.

summary(mixed_bobwhite)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: HW_Dist ~ (HuntDay * ID) + (1 | ID)
##    Data: bobwhite
## 
## REML criterion at convergence: 1969
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.85694 -0.39986 -0.09686  0.24302  3.07350 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 3679     60.65   
##  Residual             1744     41.76   
## Number of obs: 200, groups:  ID, 6
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)        16.6996    61.8382 188.0000   0.270   0.7874  
## HuntDayRabbit      30.7117    14.0785 188.0000   2.181   0.0304 *
## ID2                -3.1075    87.2444 188.0000  -0.036   0.9716  
## ID3                71.4500    87.7290 188.0000   0.814   0.4164  
## ID4               -12.6545    87.8670 188.0000  -0.144   0.8856  
## ID5                -0.6957    87.6185 188.0000  -0.008   0.9937  
## ID6               -16.6996    87.8670 188.0000  -0.190   0.8495  
## HuntDayRabbit:ID2  33.2485    19.1667 188.0000   1.735   0.0844 .
## HuntDayRabbit:ID3  -3.5830    21.1724 188.0000  -0.169   0.8658  
## HuntDayRabbit:ID4 -31.1473    23.0762 188.0000  -1.350   0.1787  
## HuntDayRabbit:ID5 -27.9975    22.5121 188.0000  -1.244   0.2152  
## HuntDayRabbit:ID6 -22.5867    22.9181 188.0000  -0.986   0.3256  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) HntDyR ID2    ID3    ID4    ID5    ID6    HDR:ID2 HDR:ID3
## HuntDayRbbt -0.167                                                          
## ID2         -0.709  0.118                                                   
## ID3         -0.705  0.118  0.500                                            
## ID4         -0.704  0.117  0.499  0.496                                     
## ID5         -0.706  0.118  0.500  0.497  0.497                              
## ID6         -0.704  0.117  0.499  0.496  0.495  0.497                       
## HntDyRb:ID2  0.123 -0.735 -0.152 -0.086 -0.086 -0.087 -0.086                
## HntDyRb:ID3  0.111 -0.665 -0.079 -0.183 -0.078 -0.078 -0.078  0.488         
## HntDyRb:ID4  0.102 -0.610 -0.072 -0.072 -0.179 -0.072 -0.072  0.448   0.406 
## HntDyRb:ID5  0.104 -0.625 -0.074 -0.074 -0.073 -0.162 -0.073  0.459   0.416 
## HntDyRb:ID6  0.103 -0.614 -0.073 -0.072 -0.072 -0.072 -0.180  0.451   0.408 
##             HDR:ID4 HDR:ID5
## HuntDayRbbt                
## ID2                        
## ID3                        
## ID4                        
## ID5                        
## ID6                        
## HntDyRb:ID2                
## HntDyRb:ID3                
## HntDyRb:ID4                
## HntDyRb:ID5  0.382         
## HntDyRb:ID6  0.375   0.384

We can check our model using the Performance package:

performance::check_model(mixed_bobwhite)

There is an unequal number of points per day for each individual covey, which makes our study unbalanced. We can check the accuracy of this model and account for this by using emmeans:

bobwhite_means <- bobwhite %>%
  group_by(HuntDay) %>%
  summarise(mean_HW_Dist=mean(HW_Dist),
            se_HW_Dist=sd(HW_Dist)/sqrt(n()))
bobwhite_means
## # A tibble: 2 × 3
##   HuntDay mean_HW_Dist se_HW_Dist
##   <chr>          <dbl>      <dbl>
## 1 Quail           22.3       4.29
## 2 Rabbit          57.0       5.30
bob_means <- emmeans(mixed_bobwhite, "HuntDay")
## Cannot use mode = "kenward-roger" because *pbkrtest* package is not installed
## NOTE: Results may be misleading due to involvement in interactions
bob_means
##  HuntDay emmean   SE  df lower.CL upper.CL
##  Quail     23.1 25.3 188   -26.92     73.1
##  Rabbit    45.1 25.1 188    -4.31     94.5
## 
## Results are averaged over the levels of: ID 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95
em_bob <- as.data.frame(bob_means)
em_bob
##   HuntDay   emmean       SE  df   lower.CL upper.CL
## 1   Quail 23.08173 25.34776 188 -26.920856 73.08431
## 2  Rabbit 45.11574 25.05572 188  -4.310749 94.54224
em_bob$HuntDay <- factor(em_bob$HuntDay, levels= c("Quail","Rabbit"))

ggplot(em_bob, aes(x=HuntDay, y=emmean)) + 
  geom_point(size=5, color="#87b8d3") +
  geom_errorbar(aes(ymin=emmean-SE, ymax=emmean+SE), width=.2, color="#6699CC") +
  geom_point(size=5, data=bobwhite_means, x=bobwhite_means$HuntDay, y=bobwhite_means$mean_HW_Dist, color = "#336699") +
  theme(axis.text.x = ggtext::element_markdown(color = "tan4", size = 12)) +
  scale_x_discrete(labels = labels) +
  theme(plot.caption=element_text(size=9, hjust=0, margin=margin(15,0,0,0)))+
  theme_bw()+
  labs(title="Comparing raw means to emmeans", caption="Light blue = raw means, Dark blue = adjusted means", x= "Hunting Day (1=Quail, 2=Rabbit)", y = "Average distance from hardwood cover (m)") 

We can see that the raw means and the adjusted emmeans don’t differ much, which lets us know that our model works for this analysis without having to account for the unbalanced points.

\(~\)

\(~\)

Nested Hierarchical Model

Dataset: Snijders, Lysanne; van Oers, Kees; Naguib, Marc (2017), Data from: Sex-specific responses to territorial intrusions in a communication network: evidence from radio-tagged great tits, Dryad Dataset

Great Tits (Parus major)

Advertisement signaling is usually linked to intersexual selection and intrasexual competition and thus is a key component of a species’ ecology. Using a novel spatial tracking system, the authors tested whether or not the spatial behavior of male and female great tits (Parus major) changes in relation to the response of a territorial male neighbor to an intruder. They tracked the spatial behavior of male and female great tits (N = 13), 1 hr before and 1 hr after simulating territory intrusions. The individual bird is a random effect, and the sex of the bird and measurements taken before and after are fixed effects.

tit <- read.csv("great_tits.csv")
tit$sex <- as.factor(tit$sex)
ggplot(tit, aes(sex, dist_m, colour = as.factor(id), shape=as.factor(measure))) + 
  geom_jitter(width =0.15, size=5, alpha=0.6)+
  ylab ("Response Distance from Playback Location (meters)") +
  xlab ("Sex") +
  annotate("text", x = 2, y = 76, label = "13 birds") +
  annotate("text", x = 2, y = 72, label = "2 measurements per bird") +
  annotate("text", x = 2, y = 67, label = "26 total measurements", size=4, color="blue")+
  labs(title="Sex-specific Responses to Territorial Intrusions", caption="1 = female
       2 = male")+
  scale_shape_discrete(name= "Measurements 
  Before / After")+
  theme_bw()

In this study, 13 individuals were measured 2 times, 1 hour before playback and 1 hour after. The individuals are separated by sex to look at the different distances between females and males. Distance on the Y axis is measured in meters, with distance being how many meters away a bird was detected before and after a playback recording was played of a territorial male’s song. Based on the measurements of distance before and after playback, it looks like there isn’t really a pattern of whether or not sex plays a role in an individual bird’s movement after a territorial playback.

Using a linear mixed model to see if there’s an interaction between the sex of the bird and distance of detection:

lmmtit <- lmer(dist_m ~ (1|sex/id), data = tit)
summary(lmmtit)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: dist_m ~ (1 | sex/id)
##    Data: tit
## 
## REML criterion at convergence: 226.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3975 -0.5381 -0.1098  0.4360  1.6636 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id:sex   (Intercept) 236.40   15.38   
##  sex      (Intercept)   0.00    0.00   
##  Residual              99.41    9.97   
## Number of obs: 28, groups:  id:sex, 14; sex, 2
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)   27.711      4.521 13.000    6.13  3.6e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
anova(lm(dist_m ~ sex/id, data = tit))
## Analysis of Variance Table
## 
## Response: dist_m
##           Df Sum Sq Mean Sq F value Pr(>F)
## sex        1   47.4   47.40  0.1309 0.7207
## sex:id     2   90.6   45.31  0.1251 0.8830
## Residuals 24 8692.3  362.18

Our p-values are not reading as statistically significant, which matches what we can see in the plot. What this can be interpreted as is that there is not a significant interaction between sex of the bird and the distance it can be detected before or after a territorial playback is used.

\(~\) \(~\)